function [ParetoFront, ParetoSet, EvaluationNum, GenerationNum] = mommop(Problem, Algorithm)
%INPUT
%	Problem.
%       Dimension: scalar, the dimension of the decision variables
%       LowerBound: 1 -by- d vector, the lower bound of the decision variables
%       UpperBound: 1 -by- d vector, the upper bound of the decision variables
%       ObjFunc: the function handle of the objective funtions,
%           y = ObjFunc(x); x: n -by- d matrix; y: n -by- 1 vector
%   Algorithm.
%       PopulationSize: scalar, the population size
%       F: the mutation factor of the differential evolution
%       CR: crossover index of the differential evolution
%       scale: the scale factor of the modified objective function
%       MaxFEs: the maximal budget
%OUTPUT
%   fval: the optimal value
%   sol: the optimal solution
%   EvaluationNum: the total number of evalutions used
%   
%%
% PROBLEM SETTING
d = Problem.Dimension;
LB = Problem.LowerBound;
UB = Problem.UpperBound;
ObjFunc = Problem.ObjFunc;
% ALGORITHM ASSIGNMENT
PopulationSize = Algorithm.PopulationSize;
F = Algorithm.F;
CR = Algorithm.CR;

N = Algorithm.MaxFEs;
NumObj = d*2;  % number of objective function
%% INITIALIZATION
GenerationNum = 1;  % the current generation number
population = rand([PopulationSize,d]).*repmat((UB-LB),[PopulationSize,1])+repmat(LB,[PopulationSize,1]);
fitness = ObjFunc(population);
EvaluationNum = PopulationSize;  % the budget used
max_fitness = max(fitness);
min_fitness = min(fitness);
population = [zeros(PopulationSize,2),fitness,modified_obj(population, fitness),population];  % [front, crowded rank, fitnesses, obj, solution]
population_merge = zeros(2*PopulationSize, 3+NumObj+d);
population_merge(1:PopulationSize,:) = population;
MutationVector = zeros(PopulationSize,d);
TrialVector = zeros(PopulationSize,d);
%% OPTIMIZATION
while EvaluationNum<N
    GenerationNum = GenerationNum+1;
    population_merge((PopulationSize+1):end,:) = 0;
    % mutation
    for i = 1:PopulationSize
        r1 = ceil(rand().*PopulationSize);
        r2 = ceil(rand().*PopulationSize);
        while r2==r1
            r2 = ceil(rand().*PopulationSize);
        end
        r3 = ceil(rand().*PopulationSize);
        while r3==r1 || r3==r2
            r3 = ceil(rand().*PopulationSize);
        end
        MutationVector(i,:) = population(r1,(4+NumObj):end)...
            +F.*(population(r2,(4+NumObj):end)-population(r3,(4+NumObj):end));
        MutationVector(i,:) = max([MutationVector(i,:);LB]);
        MutationVector(i,:) = min([MutationVector(i,:);UB]);
    end
    % crossover
    for i = 1:PopulationSize
        u1 = (rand([1,d])<=CR);
        TrialVector(i,:) = population(i,(4+NumObj):end).*(1-u1) + MutationVector(i,:).*u1;
        u2 = ceil(rand().*d);
        TrialVector(i,u2) = MutationVector(i,u2);
    end
    population_merge((PopulationSize+1):end,(4+NumObj):end) = TrialVector;
    fitness = ObjFunc(TrialVector);
    population_merge((PopulationSize+1):end,3) = fitness;
    max_fitness = max([max_fitness;fitness]);
    min_fitness = min([min_fitness;fitness]);
    EvaluationNum = EvaluationNum+PopulationSize;
    population_merge(:,4:(3+NumObj)) = modified_obj(population_merge(:,(4+NumObj):end), population_merge(:,3));
    % ranking
    population_merge = non_dominated_sort(population_merge);
    population_merge = crowded_comparison(population_merge);
    [~,rank_id] = sort(population_merge(:,2));
    population_merge = population_merge(rank_id,:);
    population = population_merge(1:PopulationSize,:);
end
ParetoFront = population(population(:,1)==1,3);
ParetoSet = population(population(:,1)==1,(4+NumObj):end);

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function y = modified_obj(x, fitness)
        scale = 40*d*((EvaluationNum/N)^3);
        xn = size(x,1);
        fitness = abs((fitness-min_fitness)./(max_fitness-min_fitness));
        fitness = repmat(fitness,[1,d]).*repmat(UB-LB,[xn,1]).*scale;
        y = zeros(xn,d*2);
        for ii = 1:d
            y(:,[2*ii-1,2*ii]) = [fitness(:,ii)+x(:,ii),fitness(:,ii)+1-x(:,ii)];
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [dominate_or_not] = dominate(fitness1, fitness2, x1, x2)
    % dominate_or_not = 1: individual1 dominate individual2
    %                   0: individual1 does not dominate individual2
        dis = sqrt(sum((x1-x2).^2));
        if dis<radius
            if fitness1(1)<fitness2(1)
                dominate_or_not = 1;
            else
                dominate_or_not = 0;
            end
        else
            obj1 = fitness1(2:end);
            obj2 = fitness2(2:end);
            dominate_or_not = (sum(obj1<=obj2)==NumObj);
            if dominate_or_not==1
                for ii = 1:d
                    if sum(obj1((2*ii-1):(2*ii))<obj2((2*ii-1):(2*ii)))==0
                        dominate_or_not = 0;
                        break
                    end
                end
            end
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [population] = non_dominated_sort(population)
    % Kalyanmoy Deb, Amrit Pratap, Sameer Agarwal, and T. Meyarivan. A Fast
    %   Elitist Multi-objective Genetic Algorithm: NSGA-II. IEEE Transactions
    %   on Evolutionary Computation, 6(2):182-197, April 2002.
    %population: [non-dominated rank, fitnesses, objective function, solution]
        population(:,1) = 0;
        population_size_f = size(population,1);
        SS = zeros(population_size_f,population_size_f);  % the solutions that be dominated by p
        NumSS = zeros(population_size_f,1);
        nn = zeros(population_size_f,1);  % the number of solutions that dominate p
        Front_f = zeros(population_size_f,population_size_f);  % the non-dominated front
        NumF = zeros(population_size_f+1,1);
        Norm_x = (population(:,(4+NumObj):end)-repmat(LB,[population_size_f,1]))./repmat(UB-LB,[population_size_f,1]);
        for p = 1:population_size_f
            for q = 1:population_size_f
                if dominate(population(p,3:(3+NumObj)), population(q,3:(3+NumObj)), Norm_x(p,:), Norm_x(q,:))
                    NumSS(p) = NumSS(p)+1;
                    SS(p,NumSS(p)) = q;
                elseif dominate(population(q,3:(3+NumObj)), population(p,3:(3+NumObj)), Norm_x(q,:), Norm_x(p,:))
                    nn(p) = nn(p)+1;
                end
            end
            if nn(p) == 0
                population(p,1) = 1;
                NumF(1) = NumF(1)+1;
                Front_f(1,NumF(1)) = p;
            end
        end
        ii = 1;
        while NumF(ii)>0
            for jj = 1:NumF(ii)
                p = Front_f(ii,jj);
                for kk = 1:NumSS(p)
                    q = SS(p,kk);
                    nn(q) = nn(q)-1;
                    if nn(q)==0
                        population(q,1) = ii+1;
                        NumF(ii+1) = NumF(ii+1)+1;
                        Front_f(ii+1,NumF(ii+1)) = q;
                    end
                end
            end
            ii = ii+1;
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [population] = crowded_comparison(population)
    % Kalyanmoy Deb, Amrit Pratap, Sameer Agarwal, and T. Meyarivan. A Fast
    %   Elitist Multi-objective Genetic Algorithm: NSGA-II. IEEE Transactions
    %   on Evolutionary Computation, 6(2):182-197, April 2002.
    %population: [non-dominated rank, crowded rank, fitnesses, solution]
        maxobj = max(population(:,4:(3+NumObj)));
        minobj = min(population(:,4:(3+NumObj)));
        max_front = max(population(:,1));
        rank_finished = 0;
        for ii = 1:max_front
            solutionId_f = find(population(:,1)==ii);
            fitness_temp = population(solutionId_f,4:(3+NumObj));
            solutionN_f = length(solutionId_f);
            distance = zeros(solutionN_f,1);
            for jj = 1:NumObj
                [~,rank_f] = sort(fitness_temp(:,jj));
                distance(rank_f(1)) = inf;
                distance(rank_f(end)) = inf;
                for kk = 2:(solutionN_f-1)
                    distance(rank_f(kk)) = distance(rank_f(kk))...
                        +(fitness_temp(rank_f(kk+1),jj)-fitness_temp(rank_f(kk-1),jj))/(maxobj(jj)-minobj(jj));
                end
            end
            [~,rank_f] = sort(distance,'descend');
            population(solutionId_f(rank_f),2) = (1:solutionN_f)'+rank_finished;
            rank_finished = rank_finished+solutionN_f;
        end
    end
end

